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Abstract 

We present a method for nonrigid registration of 2-D geometric shapes. Our 
contribution is twofold. First, we extend the classic chamfer-matching energy 
to a variational functional. Secondly, we introduce a meshless deformation 
model that can handle significant high-curvature deformations. We represent 
2-D shapes implicitly using distance transforms, and registration error is de- 
fined based on the shape contours' mutual distances. In addition, we model 
global shape deformation as an approximation blended from local deforma- 
tion fields using partition-of-unity. The global deformation field is regularized 
by penalizing inconsistencies between local fields. The representation can be 
made adaptive to shape's contour, leading to registration that is both flexi- 
ble and efficient. Finally, registration is achieved by minimizing a variational 
chamfer-energy functional combined with the consistency regularizer. We 
demonstrate the effectiveness of our method on a number of experiments. 
Keywords: meshless models, shape registration, shape correspondence, 
variational chamfer matching, distance transform 
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1. Introduction 

Registering 2-D shapes that have been deformed by nonlinear mappings is 
a challenging problem that has applications in many areas including medical 
imaging [lj and shape recognition [21 [3j. Existing shape registration meth- 
ods differ in three main components [HE]: shape representation, deformation 
model, and registration criterion. A number of recent works represent shape 
contours implicitly as zero level-sets of distance transforms [U El EU [7] . This 
shape representation has two important advantages. First, distance trans- 
form is a generic representation that can handle arbitrary shapes in arbitrary 
dimensions [5j . Secondly, distance transform provides a 2-D embedding of 1- 
D shape contours, so contours can be registered by aligning their 2-D distance 
transforms as done in image registration [H[5j[6]. Here, the shape-registration 
criterion can be simply the sum-of-squared-differences (SSD) of the shapes' 
distance transforms [4 J or mutual information (MI) [5]. Furthermore, non- 
linear shape deformations can be represented using parametric models such 
as free- form deformation (FFD) [5] and radial basis functions (RBF) [I]. 

Despite these advances, using distance transforms for shape representa- 
tion in registration methods present some common drawbacks. One issue 
is that distance transforms are redundant representations of shape contours 
which brings extra computational cost when performing registration. There 
is also the issue of stability of the representation when shapes undergo high- 
curvature deformation. Existing methods [H [5j E] use a narrow-band func- 
tion to confine computation to be around the shape contours, but it only 
partially addresses the problem and complicates the registration framework. 
An additional issue is the use of regular mesh of control points to repre- 
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sent deformations. FFD-based deformation models 0, [6] rely on a mesh of 
regularly distributed control points that makes it difficult to adapt to the 
shape contours, and the registration accuracy may suffer from folding effect 
of the control-point mesh under large deformation. The limitation of control- 
point meshes is shared by many areas including computer graphics [8j, en- 
gineering j9j [10] , and image registration [TTI [T2] . The lack of flexibility of 
regular mesh models can be addressed by replacing the mesh-based defor- 
mation models with so-called meshless models that do not rely on explicit 
connections between control points. Examples of meshless models include 
thin-plate splines using radial basis functions (RBFs) that has long been 
used for image and shape registration p21 CD, [C2, E] • However, RBFs are less 
accurate than mesh models [10], and can be computationally expensive [12j. 
A number of recent works resorted to more accurate meshless models such as 
partition-of-unity (PU) method [lOj [HI [9j [8] . The PU-meshless model rep- 
resents a continuous function by blending together local polynomial models, 
and the local models can be arbitrarily distributed and overlaid without re- 
lying on connected control points [9j. However, existing PU methods rely on 
complicated functionals (TTj [15] for regularizing shape deformation to avoid 
degenerating solutions and to produce smooth results. 

In this paper, we follow the work in [H[5j[6], an d represent shape contours 
using distance transforms. In Section [3| we briefly introduce two represen- 
tative registration methods [5j 0] using distance-transform representations, 
and elaborate on their limitations. Then, we address these limitations, and 
improve the state-of-the-art in two aspects. First, we improve the registra- 
tion criterion, by proposing a variational extension to the classic chamfer- 
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matching method [2] (Section [3]). The proposed functional does not rely on 
a narrow-band function, and can better handle high-curvature shape defor- 
mations, leading to improvements in the registration accuracy. Chamfer- 
matching has been commonly used for detecting objects under affine trans- 
formations [3j [TB] . Here, we adopt it for variational nonrigid contour 
registration. In addition, we derive a simplified gradient function of the vari- 
ational chamfer-matching functional to facilitate the registration process. 

Secondly, we adopt the PU-meshless model to represent nonlinear shape 
deformations (Section [5]). The meshless model can better represent large 
shape deformations, and can be easily adapted along shape contours to re- 
duce computational cost. Figure [T] shows the overall scheme of our method. 
To regularize the PU model, we propose a novel regularizer that penalizes 
inconsistencies between neighboring local models. The inconsistency regular- 
izer is much simpler than existing methods [HJ HE] , and provides a theoretical 
upper-bound to some regularizers [TT]. It is worth pointing out that a sim- 
ilar consistency regularizer has been used for smoothing B-splines [TT] that 
penalizes differences of adjacent B-splines' coefficients. However, this reg- 
ularizer cannot be directly applied to PU-meshless model as it ignores the 
differences of local models' coefficients caused by shifting coordinate systems. 
Our method compensates for this shifting effect using a linear operator, so 
that neighboring polynomial models can be compared on a common ground. 

In Section [6| we introduce our registration algorithm that combines the 
variational chamfer-matching gradients, the PU-meshless deformation model, 
and the consistency regularizer into a functional-minimization framework. 
We verify our method by comparing it to two recent registration meth- 
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Figure 1: Meshless shape registration, (a) Source (red curve) and target (blue curve) 
shapes. Local polynomial deformation models are defined on overlapping patches along 
the contour. Patches are weighted by locally supported radial functions. Three patches are 
shown with their corresponding weighting functions, (b) Forward and backward registra- 
tion gradients, (c) Blended global deformation map and correspondence after registration. 



ods [5j [Q. The first method by Huang et al. [5] uses distance-transform 
representation and FFD-based deformation model, while the second method 
by Chen et al. pQ represents shapes using shape-context [13] and RBF-based 
deformation model. Experiments show that our method outperforms both 
methods in terms of registration accuracy. Additionally, our method can 
handle shapes with large high-curvature deformations (Section [7]). 

In summary, our main contribution is twofold. First, we improve the 
registration criterion used by previous shape-registration methods jH El |6] 
that use distance-transform representation, and formulate the registration 
problem as the one of minimizing a variational chamfer-matching functional. 
Secondly, we introduce a meshless model for representing nonlinear shape 
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deformations, and design a simple consistency regularizer to both produce 
smooth shape deformation and avoid degenerated results. Finally, we would 
like to comment that this work extends our previous works in [THl 02] • 

1.1. Our assumptions 

We make an important assumption that shape contours are already "roughly 
aligned" before applying our registration algorithm. This is assumption is 
common to many existing shape and image registration methods that employ 
a global-to-local scheme where the shapes are registered in two stages [UEHI]. 
In these methods, shapes are first roughly aligned using affine transforma- 
tions to compensate for global shape deformations such as rotation, shifting, 
and scaling. Then, a second step seeks for local nonlinear deformations that 
align shapes as close as possible. Global shape alignment is not our main 
goal, since it can be achieved using off-the-shelf methods such as shape con- 
text [UH3], mutual information [5j, and chamfer matching [2]. As a result, 
we assume that shapes are aligned beforehand using a rigid transformation, 
and focus ourselves on the nonrigid registration part. 

2. Related work 

2.1. Shape matching and registration based on distance transform 

The distance transform (distance map) is an implicit shape representation 
whose value at a point in map indicates the minimum distance of that point 
to the shape. This implicit representation has a number of advantages. First, 
the calculation of distance transforms is highly efficient [20]. Secondly, the 
distance information helps efficient shape alignment by allowing the search- 
ing step to be adjusted according to the shapes' mutual distance [2j. Indeed, 
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chamfer-matching are frequently used for shape contour matching and de- 
tection [2T1 |3l [22j [T6] . Earlier chamfer-matching methods |2T1 [3] focus on 
detecting shapes under affine transformations with little intra-class varia- 
tions, while recent works include the use of orientation information [231 [T6] . 
hierarchical template matching [21], and statistical learning [241 1221 [25] that 
significantly improve the handling of shape variations and occlusion. The 
majority of these works apply distance transforms to shape contours, with 
exception of a recent work by Bai et al. [26] that uses distance transforms 
for matching articulated shape skeletons undergoing deformations. 

Distance transforms are also used for registering contours [23, HJ El E3, EHJ 
[29] . and point sets [30]. Similarity can be drawn between the problems of 
shape registration and shape detection, since both problems involve shape 
matching, but important differences exist. The goal of shape detection is to 
locate a given shape in an image, while registration assumes the shape's ex- 
istence and aims at accurately recovering its deformation. Paragios et al. [I] 
introduced the idea of using distance transform as an implicit embedding 
of shape contours, and registered distance maps essentially as normal 2-D 
images. Huang et al. [5j extended the method using FFD-based parametric 
deformation model. Munim et al. [7] extended the representation in [5] to 
a vectorized one for registering open shape contours. Rousson et al. [6j and 
Taron et al. [3]] investigated the integration of statistical priors into the reg- 
istration framework. In this paper, we focus on improving the registration 
criterion and parametric deformation model based on the work in [H [5] . 

Finally, distance transform is a level-set representation that can naturally 
represent shapes of arbitrary topology [5] . One may find similarities between 
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our work and level-set methods J32j|33] in three aspects. First, both methods 
use level sets of certain potential functions to represent a deforming contour. 
However, the potential functions in level-set methods are not necessarily 
distance maps, unless a regularizer is applied to the potential function to 
encourage it to maintain the shape of a distance transform [33]. In addition, 
many level-set methods aim at image segmentation (i.e., to locate a shape 
contour) , and they deform a potential function according to an external force 
that favors image features (e.g., edges). In contrast, the shape registration 
methods in [H [5j [6] and our method assume that two shape contours (i.e., 
the source and the target shapes) are given beforehand, and try to find how 
the source contour can be deformed into the target contour. Even though 
level-set methods are also used for object tracking [34j [35], their goal is to 
locate a moving object and seldom provide a dense correspondence between 
shape contours. Indeed, level-set methods can provide input data to our 
method by extracting the contours. 

A second similarity between ours and level-set methods [32], [33j I3H [35] is 
the use of penalty functions for regularization. In level-set methods [321133] . 
the regularizers penalize the geometric properties (e.g., curvature) of shape 
contours or their potential functions, to produce smooth shape boundaries. 
In contrast, registration methods [H [5] penalize the behavior of shape con- 
tours to produce smooth shape deformation, and usually do not care about 
the shapes' curvature. Finally, both level-set methods and our work min- 
imize variational functionals by iteratively solving partial differential equa- 
tions (PDEs). Please refer to [H [5] for more details of the relationship be- 
tween level-set methods and shape registration based on distance transforms. 
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2.2. Modeling and regularizing nonlinear deformations 

Distance transforms provide an implicit representation for shape contours 
themselves, but not of their deformations. It has been shown in [5] that FFD 
models can significantly improve registration accuracy over the nonparamet- 
ric deformation model [4]. The FFD model belongs to a group of so-called 
mesh models [10] that rely on a mesh of explicitly connected control points, 
that makes it difficult to adapt mesh models to different shapes [36j, es- 
pecially for shapes undergoing topology changes, where remeshing is often 
required and can be prone to errors [37]. Another drawback of mesh mod- 
els is that control points may collapse together (i.e., the folding effect) and 
change spatial relationship under large deformation, leading to numerical in- 
accuracies jTU [TO] . For mesh models, the folding effect can be addressed 
by using extra regularizers to produce diffeomorphic deformations [38], but 
these regularizers increase the complexity of the registration framework. 

To address these two issues, researchers in different areas seek for para- 
metric models, called meshless models [10] that do not rely on explicit con- 
nections between control points. Meshless models can be easily adapted to 
shape contours and greatly alleviate the folding effect of large deformations. 
Please refer to [39] for an extensive review of meshless models. Here, we 
focus on the works that are more relevant to nonrigid registration. An exam- 
ple of meshless deformation models are the Thin-Plate Splines (TPS). These 
models represent a deformation field by a linear combination of radial basis 
functions [12] centered at scattered points, without explicit connections be- 
tween them. However, it has been shown that RBFs can be numerically less 
accurate than the partition-of-unity method [10], due to the fact that RBFs 
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cannot exactly represent a polynomial function (lack of reproducibility). In 
addition, it has been acknowledged in image registration [12J and computer 
graphics [37] RBFs have higher computational cost than the PU method. 

In this paper, we adopt PU to represent nonrigid shape deformations. 
Additionally, we propose a simple regularizer to produce smooth deforma- 
tion fields and to avoid degenerated results. Existing registration methods 
regularize the deformation model using a conformity functional [H] or by im- 
posing mechanical energy constraints [ID] . There are other meshless methods 
proposed for image registration j4TJ H21 |43j [44] and segmentation [45j . Many 
of these rely on predefined landmarks and use them as boundary conditions 
to solve mechanical PDEs. Finally, there are registration methods that rep- 
resent shapes using PU models j46j[47]. These should not be confused with 
our method that use PU models to represent the shapes' deformations. 

3. Distance functions and nonrigid registration 

The goal of shape registration is to deform a source shape onto a tar- 
get shape. This is achieved by searching for the best deformation field that 
minimizes a dissimilarity measure between the shapes. Formally, if S and 
<D represent source and target shapes, respectively, and F is a dissimilar- 
ity measure between the two shapes, we seek for a warping field u(x) that 
satisfies the following equation: 

argminF(^(x),5(x / ),x / ), x' = x + u(x), (1) 
x / 

where x is a coordinate vector. The dissimilarity measure F usually depends 
on the shape model. In this paper, we implicitly represent a shape S as the 
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zero level set of its distance transform U s [U [5], where S defines a partition 
of the image domain Q. The model is given by: 



n 5 = < 



0, x G 5 

+£> 5 (x) > 0, x G R s > ( 2 ) 

-^(x)<0, x€[fi-i2 s ] 



where D 5 is the minimum Euclidean distance between location x and shape 
5, and Rs is the region inside S. Distance transforms are essentially solutions 
of a Eikonal equation given by 



||Vn 5 (x)|| = 1, x G Br, with n 5 (x) =0, x G 5, (3) 

where shape contour S defines the boundary condition. Distance transforms 
are 2-D functions themselves, and rigorously IL^x) should be written as 
II 5 ( x )(x), but to simplify the notation, we write as II 5 (x), n^( x ) and even U s . 
This simplification should not cause confusion given the context of its usage. 
The Eikonal equation in ^ is well-posed, meaning its solutions (distance 
transforms) depend continuously on the initial conditions (shape contours). 
As a result, similarity of distance transforms indicates the similarity of shape 
contours, and one can solve the shape registration problem indirectly by 
aligning shapes' distance transformations. Formally, the registration problem 
defined in ([!]) is converted to the following one given by [4]: 

argminF(II (Z) (x),n 5 (x / ),x / ), x' = x + u(x). (4) 

In this approach, registering two 1-D contours is converted to the problem 
of registering their 2-D embeddings. The later problem can be easily solved 
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using established 2-D image registration methods. In the following text, we 
will call this approach the shape- embedding method. Indeed, exactly as regis- 
tering 2-D images, F can be simply defined as the squared-sum of differences, 
and registration is achieved by minimizing following functional j5j [4] : 



E(u) = 




data term smoothness regularizer 



Here, U s and 11^ are distance transforms of the source and target shapes, 
respectively. u x and Uy are the components of the deformation field, i.e., 
u(x) = (u x (x), u y (jx)). First-order derivatives of u x and u y of large magnitude 
are penalized by a regularizer to produce a smooth deformation field and to 
avoid degenerated solutions. Finally, a proximity function N$ limits the data- 
term evaluation to be near the shape's boundary. The use of a proximity 
function is a limitation of the shape-embedding methods, and the reason will 
be clarified in Section 13.11 

Similarly to nonrigid 2-D image registration, Equation [5] has been previ- 
ously extended in several ways. First, the image deformation field u can be 
modeled parametrically using B-splines [5] , or thin-plate splines [Tj . Secondly, 
statistical priors can be leveraged to address uncertainties in the registration 
process [31J. Despite these developments, shape-embedding methods are still 
limited in some aspects. In the following section, we will elaborate on these 
limitations that motivated us for this work. 

3.1. Limitations of the shape- embedding approach 

The shape-embedding approach in Equation [5] facilitates the use of ex- 
isting nonrigid 2-D registration techniques to solve 1-D shape registration 
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problem. However, several issues need to be considered. First, as pointed 
out by El Munim and Farag the distance map defined in Equation [2] 
cannot model open shape contours or shapes that do not define closed image 
domains (i.e., "inside" and "outside" parts of a shape). This problem was 
addressed in [7J using vector distance transform. However, for the purpose of 
contour registration, we may simply use an unsigned distance transform. In 
this paper, we assume the shapes to be 0-1 encoded contours, and represented 
using unsigned distance transform as: 



where D s (x) is still the minimum Euclidean distance of position x to the 
shape contour S. In the following text, we always assume U s > 0. 



n 5 (x + u) to align with IL^x). However, deformation in the distance trans- 
form only approximates the shape deformation. Formally, n 5 (x + u) ^ 
n^x+u), where n 5 (x + u) is a direct deformation of the original shape's dis- 
tance map, while n^( x+u ) is the distance map calculated from the deformed 
shape. For example, it has been noticed by Paragios et al. [4] that scaling 
a shape is not equivalent to scaling its distance transform. Here, we would 
like to further point out that modest deformations in the shape contour (i.e., 
5(x + u)) may sometimes cause significant deformations in its distance trans- 
form (i.e., n 5 ( x+ u)), since 5(x + u) is the boundary condition for the Eikonal 
PDE (Equation [3]), its deformation will be propagated into the PDE's so- 
lution. As a result, registering distance maps may be unnecessarily compli- 
cated. Figure [2] shows the distance maps of a bending curve, i.e., II 5 ( X ) and 




(6) 



Secondly, minimizing Equation [5] essentially deforms the distance map 
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(a) Source shape n 5 ( x ) (b) Deformed shape II 5 ( X+U ) 



Figure 2: The (unsigned) distance maps of a bending curve (center in yellow). The 
distance maps are color-coded, with red and blue indicating maximum and minimum 
values, respectively. Black iso-curves are plotted to indicate points of the same distance 
values, (a) Distance map of the source shape, (b) A modest bending in the shape contour 
may cause high-curvature distortions in the distance map. 

II 5 ( X+U ), together with iso-curves indicating points of same distance values. 
The deformation of II 5 ( X+U ) has much higher curvature in certain places than 
the deformation of j(x + u) itself. As a result, when directly registering the 
distance maps, much higher bending energy is needed to cope with these high 
curvature regions, leading to registration errors or even degenerated results. 

This problem is only partially addressed in [U [5] by using a proximity 
function N$ (Equation [5]) , that limits similarity measure in the proximity 
of shape contours with distance value less than 5, where the propagated 
distortion is minimal. Intuitively, N$ defines a "narrow band" along the 
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shape contour to exclude the high-curvature part. Nevertheless, care must 
be taken to choose the narrow-band width parameter 5. On one hand, if S 
is too small, there may be little overlapping areas between the two shapes, 
and the registration algorithm can easily get trapped in local minima. On 
the other hand, if S is too large, the proximity function may fail to exclude 
distorted regions. Taron et al. [3]] proposed to use an iteratively decreasing 
5, but this approach complicates the optimization algorithm and still relies 
on this extra parameter. 

Finally, registering 2-D distance functions results into extra computation 
as the method needs to register the whole image plane instead of 1-D shape 
contours. Parametric representations [5] such as B-Splines can be used to 
improve registration efficiency, but B-Spline model cannot be easily adapted 
to the shape contours as it relies on a regular grid (mesh) of control points. 
Consequently, a lot of computation is still wasted in modeling and calculating 
deformation fields at irrelevant regions. The use of a proximity function [U 
l5l [31] further reduces the computational cost, but the formulation becomes 
more complicated. More importantly, mesh-based models may suffer from 
the folding effect when the deformation is large. Figure [3] shows a result 
obtained using the method in [5], with the deformation field represented 
using a control-point mesh of B-splines Q Large deformations in parts of 
the image cause some control points to be "squeezed" together and even flip 
their positions (folding). As a result, the deformation field no longer equals 
to the interpolation of control-point displacements, leading to less accurate 

1 Obviously, the shape deformation here can be modeled as affine global deformation, 
but we treat it as large local deformation for the sake of demonstration. 
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(a) Correspondence 



(b) Deformation field 



Figure 3: Folding of mesh model under large deformation. The images are obtained using 
method in [5]. a) Correspondence of the source shape (in blue) deformed to the target 
image (in red). The deformed source shape is in green, b) Under large deformation, the 
B-spline control points collapse together and may even flip their positions. 

or even failed registration. 

Existing registration works address the folding issue in two ways. First, 
foldings in mesh models can be removed with a stronger regularizer to pro- 
duce diffeomorphic deformation fields that are not only smooth but also 
invertible [38]. For example, the folding areas shown in Figure 2(b)| violate 



one-to-one correspondence between the original and target images (thus not 
invertible), and will be penalized by a diffeomorphic regularizer. The main 
drawback of diffeomorphic registration methods is that they are relatively 
more difficult to optimize. 

Alternatively, the folding problem can also be alleviated by replacing the 
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mesh models with so-called "meshless models" that do not rely on explicit 
connections between control points, and thus are less sensitive to the posi- 
tioning of control points. Meshless models were previously used in mechanical 
engineering [9], [10] for solving partial differential equations (PDEs), and have 
been recently introduced in computer graphics [8] and image registration [11] . 
It has been shown that meshless models can improve computational accuracy 
under large deformations [101 HI]. I n this paper, we extend a meshless model 
called partition-of-unity method [21 El HI] and adopt it for shape registration. 

In summary, the registration methods in [U [31] share two common 
drawbacks: the reliance on a narrow-band function that complicates the reg- 
istration framework, and the use of a mesh-based deformation model with 
limited adaptation to the shape contour, in addition to causing a folding 
effect under large deformations. Next, we propose a novel dissimilarity mea- 
sure based on variational formulation of the classic chamfer-matching energy 
that does not rely on a proximity function, and in Section [5| we integrate 
this energy term with a meshless parametric-deformation model that can be 
naturally adapted to the shape contours, and also reduces the folding effect. 

4. Variational chamfer-matching energy 

The limitations of shape-embedding methods are partially due to the 
fact that their distance transforms are registered as normal 2-D images. In 
this case, the dissimilarity measure takes into consideration the entire em- 
bedding 2-D space, leading to a redundancy that unnecessarily distorts the 
contour-matching quality. The narrow-band function only partially reduces 
this problem. To completely remove this redundancy, we need a dissimilarity 
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function that directly measures the shape alignment "along" the source and 
target contours without taking into consideration areas around the shapes. 
Meanwhile, we still want to use distance transform to represent the shapes 
due to its many advantages. 

We start by observing that, when the source shape S is aligned with the 
target £>, the deformed shape 5(x + u) will coincide with the zero level set 
of II^, i.e., 5(x + u) II#(x) = 0. Consequently, we may enforce alignment 
between shapes by minimizing the squared sum J Q |j(x + u) U^f rfx, which 
corresponds to the classic chamfer-matching energy function j2], used for 
affine shape registration. However, this functional can be ill-posed. For ex- 
ample, the energy function will vanish for any deformation field that shrinks 
the source shape to a single point on shape contour <D. 

We can address this problem by including a symmetric term that measures 
the distance-error between target and source shapes, similarly to the classic 
symmetric chamfer-matching energy [3J. In addition, we compensate for 
scaling by dividing the distance-error by the contours' lengths, and minimize 
the following functional given by: 

E d (u) = ^- [ |5(x + u) II 2) | 2 dx+-^ f |0(x) n 5(x+u) | 2 rfx 

A? JQ Sip JQ 

= T" / 4x + u) n|dx+^- / 2>(x) n 2 s{x+u) dx, (7) 

As Jn si<d Jn 

V v ' V v ' 

forward energy E backward energy E h 

where A s = J n 5(x + u)rfx and A® = J Q £>(x)dx are contour lengths as nor- 
malizing factors, and the second equation in ([7| holds because S is a binary 
map and 5(x) 2 = 5(x) . Additionally, E d (u) is independent of the sign 
of II, and we have assumed that 11^ > and > (unsigned distance 
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transform). The registration error is now directly measured using the shape 
contours without resorting to a proximity function as in [U El [31] . As in the 
classic chamfer matching, the use of distance transform allows for efficient 
optimization [H [5] by providing a distance-adaptive gradient function. Fig- 
ure[j](b) shows the two energy terms in our shape-registration functional. We 
now examine the minimization of our variational chamfer-matching energy. 

4.1 . Functional minimization using Euler- Lagrange equation 

A common approach to minimizing a variational functional is to calculate 
its Euler-Lagrange equation [48J. For example, if the functional is given by: 

E(u) = j L(x,u(x),u(x))dx, (8) 

then its Euler-Lagrange equation is calculated as follows: 

T dL d ( dL \ n 
J = ^:-x:U-7^ =0- (9) 



du dx \<9u(x) 

The Euler-Lagrange equation can rarely be solved directly, and it is often 
solved iteratively by converting the original problem to the equivalent PDE: 

| = J(»). do) 



In solving Equation [10[ J is used as a gradient function to facilitate iterative 
gradient-descent methods. However, calculating the gradients J from the 
variational functional in Equation [7] is not always straightforward. In these 
cases [32], we need to look for equivalent or approximating constraints to 
its exact Euler-Lagrange equation. In next section, we examine the forward 
and backward energy terms separately, and derive a simple gradient function 
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for minimizing Equation [7| We begin with an intuitive and informal deriva- 
tion, and then provide a more rigorous prove in Section 4J3 that the derived 
function can indeed minimize the variational chamfer-matching energy. 



4.2. Simplified chamfer-matching gradient 

Starting with our forward energy £^(u), its Euler-Lagrange equation can 
be written as follows: 

dV 2 d5(x + u) 
^ ~ U °« dx - ° 

where S is the source shape contour represented as a binary edge map. How- 
ever, does not exist in the classic limiting sense and is numerically 
sensitive, since edge maps S and (D are not even continuous as 2-D func- 
tions, let alone differentiate. This is the price we pay for our varitationl 
chamfer-matching formulation. Fortunately, under integration, the deriva- 
tive of 5(x + u) can be exchanged to that of according to the generalized 
derivative based on distribution theory [49j as: 
-f 



an|(x) 



<9x 



-5(x + u) rfx 



-2n 2) (x)^p5(x + u) dx. (12) 
The rationale behind this conversion is illustrated in Figure [4] by a sim- 



ple numerical example. Figure 4(a) shows a 0-1 encoded binary map of a 
shape contour, with its finite differences along x and y directions shown in 
Figure |4(b) and Figure 4(c)[ respectively. It is intuitive that these difference 



maps reproduce the negative derivative operator along the shape contour. 
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(b) ds(x)/dx 



(c) ds(x)/dy 



Figure 4: Gradients of a 0-1 edge map reproduce the negative derivative operators. When 
integrated with distance map II^, the operators are essentially transfered to H®. 

After integration with the target shape's distance map IL^x), the derivative 
operators are essentially transfered over to IL^x). In this way, we may re- 
place the original gradient operator with a "better-behaved" one on n#( x ), 
since II^( X ) is different iable. Equation 12 suggests that we may approximate 



the original Euler-Lagrange equation in (11) by a simpler constraint as: 



-2IL 



<9x 



■5(x + u) = 0. 



(13) 



Let us now examine the Equation 13 in detail. From the definition of 



distance transform, we know that 



d n g(x) 

ax 



1, i.e., 



an 



iP(x) 



ax 



is a direction 



vector that points to the increasing value of II^( X ) . Taking the negative sign 
into consideration, the left-hand side of Equation [13] is a vector pointing in the 
direction that minimizes the distance of 5(x + u) to £>(x), and the magnitude 
of this vector varies along the shape contour, proportionally to the distance 



of 5(x + u) towards the target shape, i.e., II (Z) ( x )5(x + u). Figure 5(b) shows 
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the distribution of these vectors. They resemble a force that "pulls" the 
source shape to align with the target, hence we call it the forward force. 

We now turn our attention to the backward-energy term E h . From the 
definition of E h in ([7]), we obtain: 

^=2n^ +u) ^f^(x) = 0, (14) 

where the distance map n 5 ( x+u ) is a solution of the Eikonal equation in (|3]), 
given 5(x + u) as its boundary condition, so the dependence of Il 5 ( x+U ) on 
variations in u (i.e., an ^* +u) ) is difficult to analyze. This dependence can be 
greatly simplified if we assume that the deformation field is a simple shifting, 
i.e., u(x) = Uo, Vx G f2 being a constant function. With this assumption, 
the shape deformation equals to the deformation of its distance map jl]: 



n^(x+u) = n^( x + u). Substituting this into (14), we obtain: 

dL h <9II 5 (x + u) . . 

= 2 n, x+u) ^±^ (x) 



= 2II 5(x+u) — gp©(x). (15) 

We would like to clarify that ^ n ^ +u) in (15) is simply the gradient of II 5 ( X+U ) 
without further dependence on the deformed shape 5(x + u). Using a more 
complex notation, this gradient term should be written as dU ^^ where S f = 
s(y + u) is the deformed shape. Here, we keep its notation simpler and 
consistent with Equation [l3j It is also worth pointing out that the above 
derivation is valid only under the assumption that u(x) is constant for every 



x, but we relax this assumption to arbitrary continuous u(x), and use (15) 



as the gradient function for minimizing the backward energy. This is of 
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course only an approximated function as we have ignored variations of the 
deformation function u(x), but we show in the following section that this 
simplification is not arbitrary and the approximated gradient function will 
lead to the minimization of our chamfer-matching functional. 



The gradient function in Equation [15] has similar intuitive motivation 

we 



as the forward force. By negating the left-hand side of Equation 15 



obtain — ^ that acts like a force "pulling" the target shape towards the 



source as shown in Figure 5(c)[ In actual minimization, the target shape is 



kept unchanged (static), thus the force is applied in reversed direction to the 
source shape, and we call it the backward force. 



Combining (13) with (15) and ignoring the constants, we may now ap- 
proximate the Euler-Lagrange equation of the original variational functional 
(Equation [7| as follows: 

j = -n^ (x) v 5(x + u) + n 5(x+u) — ^ — ^(x). 16) 

ax ax 



forward force backward force 

J closely resemble the gradient of classic chamfer-matching method [2], and 
we call J the variational chamfer-matching gradient. J can also be inter- 
pret ated as a force field deforming the source contour, in a similar analogy 
used in the deformation of active contours [32J. Finally, the strength of the 
combined force increases according to their mutual distances. 

Note that, in the chamfer-matching energy functional (Equation [7]), we 
may also use the L 1 norm instead of the squared-sum (i.e., L 2 norm). In this 
case, the forward and backward forces degenerate to direction vectors with 
uniform magnitude regardless to the mutual distance between the shapes. 
Our experiments showed that L 1 norm is more sensitive to local minima, 
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(a) (b) (c) 

Figure 5: Variational chamfer-matching gradient, a) Source and target shape contours, 
b) Forward force is distributed along the source shape (curve in yellow). The distance 
transforms are color-coded, with red and blue indicating the maximum and minimum 
values respectively. The forward force pulls the source towards the target, c) Backward 
force in reversed direction that pulls the target to the source. The magnitude of these 
forces are proportional to the distances between two shape contours. 

and leads to slower minimization convergence. This observation echoes a 
similar finding in the classic chamfer-matching method [2]. 

4-3. Variational chamfer-matching based on the simplified gradient 



We have derived an approximate but simple gradient function J in (16) 
to minimize the variational chamfer-matching functional £^(11). Now, we 
show that J can indeed lead to the minimization of E d (u). From variational 
calculus, we know that a necessary condition for minimizing E d (u) is its 
Euler-Lagrange equation J = 0. However, this condition is not sufficient as 
the Euler-Lagrange equation may also hold when E d (u) is at local minima. 
Nevertheless, the Euler-Lagrange equation is widely used in computer vision 
to minimize variational functionals. Here, we wish to establish a similar 
relationship between the minimization of E d (u) and our approximate gradent 
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function J, then the use of J is theoretically justifiable. It turns out that 
J = is indeed a necessary and sufficient condition for the variational- 
matching energy to vanish, i.e., E d (u) = 0. 

Before elaborating on the proof, we want to clarify that our discussion 
is limited to continuous deformation fields, i.e., u G C(Q) with C(Q) as the 
set of continuous functions defined on domain f2. Vanishing of the chamfer- 
matching energy, i.e., min E d (u) = can only be reached when two shapes 

ueC(Q) 

5(x + u) and £>(x) can be exactly aligned, meaning that there exists a contin- 
uous deformation u(x), such that 5(x + u) = £>(x), Vx G A mismatch of 
the two shapes at any image pixel p G Q would lead to E d (u) > 0, due to the 
non-negativeness of unsigned distance transform. As a result, if the source 
and target shapes have different topologies, then E d (u) > 0, Vu(x) G C(Q) 
since there is no continuous mapping between shapes of different topologies. 
By focusing on E d (u) = 0, we have excluded the cases of aligning shapes 
with different topologies. Fortunately, these cases are not very interesting 
for many practical applications. When the two shapes have the same topol- 
ogy, then £" d (u) = is indeed the global minimum of E d (u). In summary, 
we can prove the following lemma: 

Lemma 4.1. // min E d (u) 0, then E d (u) = J(u) 0. 

ueC(Q) 

Proof On one hand, E d (u) = means that a continuous deformation field 
u(x) exists such that 5(x+u) = £>(x), and II^( X ) = II 5 ( X+U ). By the definition 
of distance transform, we obtain II 5 ( X+U ) £>(x) = 0, ILp( x )5(x + u) = and 

2 Rigorously, exceptions may exist for singular point set of zero measure. As far as 
discrete image data is concerned, this assertion is valid. 
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J(x) = 0. Consequently, J(x) = is indeed necessary for E d (u) = 0. 

On the other hand, let us assume J(x) = 0, Vx G fi, then we can show 
that j(x + u) = £>(x) must be true. For example, if j(x + u) ^ £>(x) for 
some x G fi, we can assume £>(x) = 1 and 5(x + u) = without loss of 
generality. As a result, the forward force vanishes at x, from the definition of 



J(x) (Equation 16). Since J(x) = 0, the backward force also become zero at 



x, i.e., — II^x) dIi Q^ ) — 0, but we have assumed that £>(x) = 1, so ILp(x) > 
1, leading to a contradiction. Consequently, J(x) =0 ensures 



d n g(x) 



and ,, x 

that 5(x + u) = £>(x) holds for all x G £1 This in turn proves E d (u) = 0, 
and hence the sufficiency of J = for E d (u) = □. 

We have shown that J(u) can be used to minimize the variational chamfer- 
matching energy. Now, we can combine the chamfer-matching energy with 
different regularizers and deformation models for nonrigid shape registration 
such as the derivative-based regularizer of Equation [5] and B-Spline repre- 
sentation in [5j. However, B-Spline models rely on an explicitly connected 
control-point grid (i.e., mesh), making it difficult to adapt the model to 
shapes of challenging topologies (holes, splits). Instead, the whole image 
domain is taken into consideration even though a large part of it can be 
irrelevant. Additionally, B-Spline models may suffer from the folding ef- 
fect that control points collapse together under large deformation, causing 
numerical instability. Next, we study how to remove the redundancy from 
parametric deformation models, by adopting a meshless representation that 
approximates the shape's deformation field. 
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5. Meshless deformation model 

Our model is derived from recent developments in computer graphics [8] 
and mechanical engineering [10J, on building shape functions of arbitrary 
topology. Specifically, our model belongs to a group of so-called partition- of- 
unity meshless methods. Although there are other meshless shape-deformation 
models based on thin-plate splines and radial basis functions (RBFs) PQ, 
those methods are less accurate than polynomial-based representations as ra- 
dial basis functions cannot exactly represent polynomial deformations (lack 
of reproducibility) [10J. In the partition-of-unity model, one can ensure poly- 
nomial reproducibility by locally modeling shape deformation as polynomi- 
als. These local deformation models are then blended together into a global 
deformation field based on a set of weighting functions. We illustrate the 
partition-of-unity blending concept using the 1-D example shown in Figure [5| 
The figure shows 1-D patches represented by their weighting functions parti- 
tioning the data domain into three regions. The local deformation models are 
represented by the colors pink, blue, and red. The global deformation model 
is represented by the green polynomial curve obtained by blending the local 
models. Here, we would like to remind the readers that the meshless model 
in our method is used for representing the incremental nonrigid deformation 
(i.e., function u in Equation [7]) , but not the shape contours themselves (the 
shapes are represented by distance functions and binary maps). Our ap- 
proach largely follows the partition-of-unity used in computer graphics [8]. 
In the following subsections, we first introduce the local deformation model, 
and then explain how to blend them into a global model. 
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local polynomial models 
blended global model 




meshless model 



1 r 




weighting functions 



Figure 6: A 1-D meshless model based on partition-of-unity. The local polynomial models 
are blended together (top), using predefined weighting functions (bottom). 



5.1. Local polynomial model 

We commence by modeling shape deformation locally in scattered local 
domains called patches [9] . As in [8] , we define these patches as the support 
of some weighting functions. Various types of weighting functions exist [10J. 
In this work, we use the 2nd-order B-Spline w p (x) = a p Wb ^ x ~ p ^ for its 
computational efficiency, where a p G (0, 1] is the patch's predefined influence 
factor in the final global blending. Formally, a 2nd-order B-Spline function 
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is defined as a piecewise polynomial: 



w(r) 



| - r 2 , r < \ 

3 



,r > § 



This weighting function's support is a disk with radius r p centered at p (ball 
in 3-D0). To summarize, a disk-like patch p is defined by three parameters 
(p, r p , i.e., its center position, radius, and influence factor. The union of 
these patches covers the image region of interest, i.e., Q C U^supp (w p ), 
p e P, with P is the set of all patches. Note that although the weighting 



function in (17) is a radial function, its usage is different from previous RBF 
models such as thin-plate splines pQ. Here, RBFs are used for blending local 
polynomial models instead of directly representing the shape deformation. 

It is also important to point out that the patches simply define a partition 
of the computational domain, and remain static during shape deformation, 
unlike B-spline control points, whose displacements directly control the de- 
formation field. Instead, the local deformation u p = (u, v) within each patch 
is defined by the parameters of a polynomial model centered at p as follows: 

s,t=m s,t=m 

u(x) = a s,tx s y l and v(x) = ^ b s , t x s y\ (18) 



Alternatively, u p (x) is a linear combination of monomial basis functions 
T (x) = (1, x, y, xy, x 2 , y 2 ) . . . , x m y m ) ) with a two-column coefficient matrix 



3 A 3-D extension is straightforward. 
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dp as: 

u p (x) = d p T p (x), (19) 

where vec (dp) = (a , , &o,o, ■ ■ ■ ,a m5m ,6 m5m ) T and p (x) = 0(x - p) is the 
basis vector re-centered at p. The sequence of monomials in <fi can be rep- 
resented as a one-to-one mapping from the natural number to nonnegative 
2-D monomial indices ^ : N ^ No x No. In this paper, the monomials are 
arranged in a Pascal-triangle manner [TO] , where low-order monomials are 
arranged before high-order ones. For monomials of the same order, the ones 
with more balanced indices in x and y are arranged first. This special ar- 
rangement increases numerical stability according to existing literature [10J. 

5.2. Blending local models into a global deformation field 

Once the local deformation models are at hand, the deformation at any 
point x is obtained by blending local fields of patches that contain x. These 
patches are denoted by AT X = {p | x G suppf^)}. The blended global- 
deformation field is given by: 

u(x) = ^( X H( X )> with r p (x) = Wp ^ . (20) 



is es- 



Here, r p (x) ensures the partition-of-unity (PU), so that Equation 20 
sentially a weighted average. Compared to the B-spline model, the partition- 
of-unity method does not need explicit connections between neighboring 
patches, and allows for arbitrary overlapping of local models. This blending 
scheme has been previously used in the polyaffine model for diffeomorphic 
image registration [50], meshless image registration [llj, and point-based 
computer graphics [8j . Some previous works [HI [9j HO] directly define the 
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partition functions r p (x) without constructing them from weighting func- 
tions. However, weighting functions provide two appealing features. 

First, weighting functions are intuitive to handle, they allow us to assign 
different weights and scales to patches, and can potentially be made adap- 
tive according to prior information of shape deformatioi^J It is desirable that 
the positions, scales, and weighting factors of patches can be optimized ac- 
cording to registration results in a data-driven fashion. However, data-driven 
approaches are often application dependent, so this topic is beyond the scope 
of our paper. Here, we exploit a simple heuristic to register shape contours. 
Broadly speaking, we choose the patches' scales according to the shapes' 
mutual distance, since edge points that are far apart are more likely to have 
large-scale deformations. We found that this simple heuristic works well in 
most cases. Figure |7(a)| shows an example shape with patches distributed 
along its contour. Figure [7(b)] shows the sum of their weighting functions (i.e., 
the denominator of partition functions r v in Equation 20). Thus, we limit the 



computation of the deformation model to take place along shape contours, 
and achieve a similar effect as the narrow-band function used in H] . 

Secondly, weighting functions also provide a way to specify our confi- 
dence about the approximation accuracy of different local deformation mod- 
els. This confidence measure will be used to regularize the global deformation 
field. Regularizes are important for preventing degenerating solutions, since 
the solutions of shape registration may not be unique. Regularizes remove 
ambiguous solutions by penalizing undesired fluctuations in the estimated de- 



4 In the partition-of-unity literature [9], other a priori knowledge can also be included 
about local behavior of the solution in the finite element space. 
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(a) 



(b) 



Figure 7: 2-D meshless deformation model for shape registration, a) Circular patches are 
distributed along the shape contours with centers marked as black crosses, and their scales 
are adapted according to distance of the source (in red) to the target shape (in blue), b) 
The sum of weighting functions, whose support covers the interested image region. 

formation field, and prefer "smooth" solutions. Next, we study the problem 
of how to regularize a meshless shape-deformation model. 

5.3. Smoothness regularizer 

We have shown that a global deformation field relating two shapes can 
be obtained by blending local deformation fields using Equation [20| A clas- 
sic regularizer used in B-spline models [5j is to penalize the magnitude of 
derivatives of the resulting deformation. However, calculating derivatives of 
our meshless model is not straightforward. For example, according to the 
product rule, first-order derivatives of the global deformation are obtained 
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as follows: 



Vu(x) = (Vr p (x)u p (x) + r p (x)Vu p (x)) . (21) 

The above differentiation involves derivatives of local deformation models 
u p (x), and those of the partition functions r p (x) given by: 

_ , v _ V^( x ) Ep/e^x ^( x ) ~ ^p( x ) E p /g^x ( x ) m ^ 
Vr p (xj — 2 . \ZZ) 

It is clear that calculation of Vu(x) is complicated by the presence of blending 
functions. This problem was partially addressed in [11], by approximating 
Vu(x) using the differentiations of local models (i.e., Vu p (x)) as follows: 

Vu(x) = ^(x)Vu p (x), (23) 

that is basically a weighted average of the derivatives of neighboring patches. 
The approximated derivative Vu(x) is then penalized in a similar fashion 
as in the classic smoothness regularizer [5J. However, this approximation 
is only accurate when neighboring patches have very similar derivatives (in 
conformity). This condition is achieved by penalizing a Sobolev conformity 
term between neighboring local deformation models as follows [TT] : 

S r = E / ™ P (xH(x)||DX(x) - DX(x)||^x, (24) 

where D ,? u p (x) = d^d^u p (x), r) x ,r) y = 0,1,..., k, and 77 = (rj x ,rj y ) is a 
multi- index vector with total order \rj\ = rj x + i] y . The usage of multi-index 77 
is to include various orders of derivatives. In other words, the Sobolev confor- 
mity term is essentially a weighted squared sum of difference (SSD) between 
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the derivatives of neighboring polynomial models. By enforcing S^ q — >■ for 
each pair of patches (p, g), the local deformation models are forced to have 



similar derivatives. Equation [24] can be simplified into a quadratic term of 
local polynomial coefficients, but this approach is still much more involved 
than regularizing mesh models, and it involves two components: the Sobolev 
conformity term and a classic gradient-based smoothness term. 

We found that a much simpler regularizer can be derived for smoothing 
the deformation field. First, the conformity between local models can be 
simply enforced by shifting and comparing the coefficients of local polynomial 
representations, instead of using an integral SSD form as in Equation [24} In 
this way, the conformity regularizer is greatly simplified, without having to 
compare the derivatives of local models. To differentiate with the conformity 
regularizer, we call our simplified functional as consistency regularize^ as 
it penalizes the inconsistencies in the polynomial coefficients of neighboring 
patches, instead of measuring the conformity of their derivatives. Secondly, 
our consistency term also penalizes fluctuations in the deformation field, 
and is enough to ensure its smoothness, further simplifying the registration 
framework by replacing the classic derivative-based regularizer. 

In Figure [8j we show the idea of our consistency regularizer using a sim- 
ple 1-D example of approximating scattered data points. Intuitively, the 
fluctuations in the global model is related to the inconsistencies among lo- 
cal deformation models (Figure |8(a)| ). One may consider that consistency 
between two local deformation fields, u p and u g , can be measured as the 
Euclidean distance between their parameters d p and d q as follows: 

E m = ||vec(d p -d g )|| 2 . (25) 
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Figure 8: Consistency regularizer for 1-D scatter-point approximation. Data points were 
generated from sinusoid function with added noise. Three patches of equal radius and 
weighting factors were used to partition the computational domain [— tt, it]. The parameter 
A controls the relative contribution of the regularizer term, a) Fluctuations in the global 
model are related to inconsistencies between local polynomial models, b-d) By increasingly 
penalizing inconsistencies among local models, the global model becomes smoother. 



Here, we stretched the parameter matrices d p , d q into vectors vec(d p — d q ). In 
fact, a similar idea has been used for smoothing B-splines [T7j by penalizing 
differences of adjacent B-splines' parameters. However, for meshless models, 



the consistency regularizer based on (25) incorrectly compares local defor- 



mation coefficients from different local coordinate systems. We demonstrate 
this effect using another 1-D example as shown in Figure |9) In this example, 
three linear (i.e., first-order polynomial) deformation models are blended into 
a global one. Even though the local models are themselves consistent, their 
coefficients are different due to different coordinate systems. The regularizer 
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local polynomial models 
blended global model 




h{x) = 0.3x - 0.5 f 2 {x) = 0.3x f 3 (x) = 0.3x + 0.5 
dj = (-0.5, 0.3) dj = (0.0, 0.3) dj = (0.5, 0.3) 



Figure 9: A 1-D example of shifting local coordinate systems. Three local models are 
blended into a global one. The local models are represented using first-order polynomials, 
with basis functions </>(x) = (l,x) T . Shifting the coordinate system leads to a linear 
transform of local model coefficients. 



in Equation [25] failed to take the coordinate transform into consideration. To 
address this issue, we need to shift and align the local coordinate systems so 
that d p and d q is compared on a common ground. 

Fortunately, aligning the local coordinate systems can be achieved using a 
linear operator. To begin with, the representation of u g (x) in the coordinate 
system of another patch p can be obtained by shifting u g (x) as follows: 

u 9 (x) = d/0 g (x) 

= d/0(x - q) 

= d/0(x - p + p - q) 

= d/0 p (x + p-_q) (26) 

Ax 
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As a result, the alignment of local deformation models reduces itself to the 
shifting of the basis functions </> p (x). Our choice of monomial basis functions 
makes this computation particularly simple. In general, shifting a polynomial 
basis 0(x) by Ax = (5x, Sy) leads to: 

0(x + Ax) = (1, x + Sx, y + Sy, (x + Sx)(y + Sy), . . . , (y + 5y) m ) T 

= S T (Ax)0(x), (27) 

where S T (Ax) is the linear basis-shifting- operator. For example, for second- 
order basis vector </>(x) = (1, x, y, xy, x 2 , y 2 ) T , its shifting-operator is: 

\ 



S T (Ax) 



V 
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Sx 
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SxSy 


S 2 x 
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(28) 



/ 



By substituting (27) into (26), we obtain the shifted local deformation field: 



iWx) = d/S T (Ax)0 p (x) 
= (S(Ax)d g ) T p (x). 



(29) 



Here, we used u g ^ p (x) to distinguish u q from its shifted representation. 
Therefore, shifting the local coordinate system leads to linearly transformed 
polynomial coefficients d q ^ p = S(Ax)d g . Now, we can modify the consis- 



tency term in Equation 25 to include the shifting operator as: 



|vec(d p - S(p - q)d q 



(30) 



37 



Finally, we can sum up the consistency term for every pair of neighboring 
patches as follows: 

EC = ^Y,(Y. w «(\\p-^ E pX (31) 

p \ q / 

where TV is the number of patches, and the terms E° are weighted by w q 
based on the distance between patches. The inclusion of weighting func- 
tions has two purposes. First, it effectively limits the interaction range of a 
patch p to its neighboring patches q that contain the center of p, otherwise 
^g(l|P — q| | ) = . Secondly, the weighting function gives more emphasis to 
the consistency of neighboring patches that are close to each other. 

5.4- Verification of the regularizer 

We have introduced a simple consistency regularizer based on the linear 



shifting operator in (31). Before adopting it for shape registration, we would 
like to briefly verify its ability of smoothing the deformation fields. We have 
shown in Figure [8] an example of 1-D approximation of scattered data points, 
combining the consistency regularizer with a SSD data-fitting cost function. 
Formally, the combined data-fitting functional is given by: 

i 

where (x^,/^(x^)) are scattered data points, and A is the weight of the con- 
sistency regularizer. Figure [8] illustrates how the increase of the regularizer 's 
weight leads to smoother global approximation, and the global model ap- 
proximates a polynomial fitting asymptotically. The effect of the consistency 
regularizer can also be analyzed theoretically. Here, we aim at providing an 
intuitive understanding of the regularizer rather than a rigorous proof. 
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On one hand, we can show that minimizing the consistency term in (25) 
can also enforce the conformity between local deformation models, i.e., the 
continuity in their derivatives. To verify, note that the differential operator on 
polynomials can again be represented as a linear operator on the coefficients 
dp. Let us denote a general differential operator and its matrix representation 
using the same symbol D v . Combining it with our shifting operator into the 



conformity term in Equation 24, we obtain: 



S r = E / ™ P (xK(x) ( D X(*) - D%^(x)) 2 dx 

= X) / ™p( x K( x ) ( D " d ^( x ) - D " S (P - q)d,0,(x)) 2 rfx 

\ v \<k J 

< ||vec(d p -S(p-q)d g )|| 2 ^ ||D1 2 /^ p (x)^(x)0j(x)dx. (33) 

Consequently, the conformity term is upper-bounded by our consistency reg- 
ularize^ and minimizing the consistency regularizer will similarly enforce the 
conformity condition among local models. 

On the other hand, when E c 0, the blended global model is coerced into 
a polynomial one, and the order of this polynomial is the same as the basis 
functions. For instance, if we choose the linear basis vector = (l,x,?/) T , 
then E c = implies that u(x) is a piecewise linear function and ||Vu x || 2 + 
|| Vu y || 2 = 0. As a result, our consistency regularizer has the same asymptotic 
property as the classic gradient-based regularizers. 

We have shown the relationship of our consistency regularizer with both 
the conformity term used in a recent partition-of-unity method [TT] and the 
gradient-based regularizers for mesh models [5]. Next, we formulate and 
solve the nonrigid shape registration problem, by combining the meshless 
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deformation model with our variational chamfer-matching energy. 



6. Registration using gradient-descent method 

With the meshless deformation model, the registration problem is now 
converted to the one of seeking the "best" deformation parameters d p) p G P, 
by combining both the variational chamfer-matching and the consistency reg- 
ularizer. The chamfer-matching energy measures the closeness of registered 
shape contours (data term), and the consistency regularizer penalizes fluc- 
tuations in the deformation field (smoothness term) . Formally, we formulate 
nonrigid shape registration as the following functional minimization problem: 

d p = arg min E° = arg min (E d (u) + \E C ) , (34) 

dp dp 

where parameter A defines the relative importance of the consistency term. 



As in previous works j5j [2] , minimizing ( 34 ) can be achieved using gradient 

descent. The gradients can be derived as follows given by: 

dE v _ dE d (u) x dE c 
dd p dd p dd p ' 



Gradients of the consistency energy can be calculated from (31) as: 

= jf E w «(IIp - ill) ( d p - s (p - iH) • ( 36 ) 

q 

Gradients of the variational chamfer-matching energy can be calculated from 

chain-rule of variational calculus, to obtain: 

dE d (u) f 9u(x) 



/ J( x) 



ddp J ddp 



-dx 



J J(x) • r p (x)0(x)rfx. (37) 
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Algorithm 1: Shape registration based on gradient-descent. 
Input: Binary maps for source and target shapes (s and £>), 

partition-of-unity patches P, consistency weighting factor A. 

Output: Deformation field u(x). 

Initialize local deformation coefficients d° <— 0, for p G P; 

Calculate distance map of the target shape II^; 

while Not converge do 

Calculate global deformation u k (x) from local deformation 

coefficients (Equation 



20); 



Warp source shape as 5 fc (x) = 5(x + u h ) ; 
Calculate distance map II 5 fc( x ); 
Update Jacobian matrix (Equation 



36 



and Equation 



37); 



Update dp +1 using gradient-descent rules; 



end 



Here, J(x) is left-hand-side of the Euler-Lagrange equation of i? d (u), and we 
use its approximation J(x) defined in Equation 



16 



Using the derived gradients, we implemented an optimization algorithm 
based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [5T| . At 
each iteration of the algorithm, the source shape is first warped using the 



deformation field reconstructed from local field parameters (Equation 20), 
and then its distance transform II 5 ( x+ll ) is updated. Both the destination 
shape <D and its distance transform 11^ remain constant. Conceptually, the 
gradient-descent process is described in Algorithm [T] 

Direct warping of a shape contour might cause deterioration of image 
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quality since the deformation 5(x+u) is often implemented using interpolated 
backward mapping. Evidently, a warped 0-1 edge map may produce decimal 
pixel values, i.e., no longer a binary image. In practice, the deterioration of 
edge maps is rarely an issue for our algorithm, and can be naturally handled 



by our variational chamfer-matching gradient (Equation 16), by relaxing the 
binary edge map to a gray-scale one. Furthermore, we can use an extended 
version of the distance transform defined on gray-scale images that takes the 
gray-scale edge map value into account. Nevertheless, in this paper, we want 
to keep the edge map binary, for better evaluation of the registration results, 
since the compared methods in [5[ d] are based on binary representations. 

In previous implementation of distance-based registration [5j [^J the shape 
contours are represented as connected line segments, that is basically a vec- 
torized representation. Warping of the shape contours is achieved by deform- 
ing the end points of these line segments. Thus, interpolation and resulting 
deterioration of edge maps can be avoided. This approach may be adopted for 
our method, but is not straightforward to implement. Fortunately, in many 
applications, the input shape contours are detected from grayscale source 
images. In these cases, we simply warp the source images, and then apply 
an edge-detection step to recover the deformed binary contour. To summa- 
rize, we modify the registration Algorithm [T] to obtain an edge-preserving 
one given by Algorithm [2] in which the main differences are highlighted in 
italic. Both the distance transform and edge detection are highly efficient 
operations, with minimal influence on the algorithm's computational cost. 



5 Available for download at http://www.cse.lehigh.edu/~huang/downloads.htm 
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Algorithm 2: Gradient-descent with edge-preserving. 
Input: Source and target images (I s and I t ) ) partition-of-unity 

patches P, consistency weighting factor A. 

Output: Deformation field u(x). 

Initialize local deformation coefficients d° <— 0, for p G P; 

Detect target shape contour from the image <D <— I t ; 

Calculate distance map of the target shape II^; 

while Not converge do 

Calculate global deformation u fc (x) from local deformation 

coefficients (Equation 20); 



Warp source image as I* (x) = 7 s (x + u 
Detect source shape contour S k <— 1% ; 
Calculate distance map H s k^\] 
Update Jacobian matrix (Equation 



36 



and Equation 



37); 



Update d^ +1 using gradient-descent rules; 
end 



Finally, we adopt the hierarchical multi-scale registration strategy used in 
previous nonrigid registration methods (5j [SU E] in order to avoid local min- 
ima (i.e., a coarse-to-fine approach). Briefly, multi-scale registration starts 
with sub-sampled images (at the coarser level). The deformation obtained 
at coarse scale is then used to initialize the algorithm at finer scales. 

7. Evaluation 

In this section, we evaluate our method on a number of experiments, and 
compare it with the methods by Huang et al. [5] and Chen and Bhanu [Tj. As 
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described in Section [3| Huang's method is based on distance-map represen- 
tation and a B-spline deformation model, while Chen and Bhanu's method 
is based on shape-context and a thin-plate deformation model. Our goal 
in this comparison is to demonstrate two aspects of our method. First, by 
removing the narrow-band function, our chamfer-matching functional can 
help improving the registration accuracy of distance-based methods at high- 
curvature regions. Secondly, by replacing spline-based mesh models with our 
partition-of-unity model, we can handle challenging shape deformations that 
may be difficult for existing methods. We tested our method on the Brown 
University shape dataset [52j, that was also used by many previous meth- 
ods j5j [U H|. Additionally, we demonstrate an application of our method to 
image sequences of cell mitosis. 

7.1. Brown University shape dataset 

Our method requires some parameters to be set. First, we empirically 
set the consistency regularized weighting factor A = 0.001. Secondly, the 
representation power of our meshless model depends on both the density and 
distribution the image domain partitions. To simplify, we began with regu- 
larly distributed patches with radius r = 20 (pixels) and the spaces between 
patches d = 6. We normalized the images of Brown University shape dataset 
to sizes 150 x 150 pixels. Finally, we globally aligned the shapes using the 



rigid registration method implemented in [53]. Figure 10 shows registration 
results produced using our method. For a quantitative evaluation, we used 
the average mutual distance between the registered shape contours [lj. For 
qualitative comparison, we also visualize the results in a similar manner as 
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(a) (b) (c) (d) (e) 



Figure 10: Brown university shape dataset. (a) Target images, (b) Overlaid target (in 
black) and source images (in red) before registration, (c) After registration, (d) Corre- 
spondence between target and source images, (e) Deformation fields as distorted grids. 
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in HI [I]. As in [lj, we selected three different shapes (persoi^J hand, and 
fish). The average pixel distances after local registration for person, fish, 
and hand were 0.14, 0.24, and 0.08, respectively. This result was better than 
the one reported in pQ (0.59, 0.45 and 0.61, respectively). Additionally, for 
most cases, the maximum pixel distance was less than 2 pixels showing that 
registration quality was consistent along contours. We believe that the use 
of both the chamfer-matching functional and meshless model contributed 
to the improved accuracy. On one hand, chamfer-matching functional is a 
more direct measurement of shape contour's alignment than shape-context 
descriptors. On the other hand, the partition-of-unity model is numerically 
more accurate than thin-plate splines that use radial-basis functions. 

7.1.1. Shapes with occlusion and large deformation 

Our method was able to register shapes with occlusion, large deformation, 
and high curvature. Here, we compared our method with Huang's method [5]. 



Figure 11 shows the registration of occluded hand shapes. In [5j, landmarks 
were placed at occluded parts to guide the registration process. Without 
landmarks, it is difficult for splines to handle large deformation from the 
source (blue) to the target (red) shape due to foldings in the mesh model 



(Figure 11(e)). In contrast, our meshless model alleviated the folding effect 
(Figure |ll(f)[ ). Here, we plotted the deformation field by deforming an image 
of grid, but it does not mean that our model is restricted to such a topology. 
Quantitatively, our method still achieved remarkable accuracy with the mean, 
maximum, and variance of edge distances as 0.15, 1.41 and 0.13, respectively, 



6 Named dude in the original dataset. 
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(a) Huang et al. (b) Ours 




(c) Before (d) After (e) Huang et al. (f) Our method 



Figure 11: Occluded shapes (hand). The contours are source shape (blue), target shape 
(red), and deformed source shape (green), a) Huang et al. [5], average distance: 0.64, max: 
9.2, variance: 1.99 . b) Our method, average distance: 0.15, max: 1.41, variance: 0.13. 
c) Deformed source (green) and target (red) shapes, d) The deformation field represented 
using a deformed grid image. 

in comparison with Huang's result (mean: 0.64, max: 9.2, variance: 1.99). 
The mesh model's folding effect may be further complicated by the narrow- 



band function N$. Figure 12 shows registration of occluded shapes (bunny) 
using increasing narrow-band S. On one hand, we observed that increasing 5 
will generally reduce the ability of Huang's method to handle high-curvature 



regions (Figure [l2(a)1 and [12(b)] ) as using large S may include embedding 
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(a) 5 = 10 



(b) 6 = 50 



(c) Ours 



Figure 12: Occluded shapes (rabbit), a) Huang et al. [5 j, average distance: 0.62, maxi- 
mum^. 38, variance: 0.89. b) Increasing the narrow-band width futher reduces the reg- 
istration accuracy, with average distance: 1.00, maximum: 10.44, variance: 2.63. c) Our 
method, average distance: 0.41, maximum: 2.23, variance: 0.28. 



space that undergoes high-curvature bending (see Figure [2] for demonstra- 
tion), while our method was not affected at all by such a parameter (Fig- 



ure (12(c)). On the other hand, with decreasing 5, Huang's method became 



prone to local minima. Figure [13] shows an example of shapes with both 
high curvature and large deformations. Here, we tested Huang's method us- 
ing different narrow-band width S. When S was small (Figure [13(a) ), the 



algorithm converged prematurely, unable to extend the deformation large 
enough. When S was large (Figure [l3(b)D , the method was unable to cope 
with high-curvature regions around the bending arm. An optimal result was 
obtained around 5 = 20. However, even this optimal result was still less 
accurate when compared to our method. 
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(a) 6 = 10 



(b) 6 = 50 



(c) 6 = 20 



(d) Our method 




Figure 13: Large deformation (dude), a-c) Method of Huang et al., with varying narrow- 
band width. Best result was obtained when 6 = 20, with average distance: 0.17, max: 
3.61, variance: 0.20. d) Ours, average: 0.14, max: 1.00, variance: 0.12. 

7.2. Adapted partition- of -unity 

We now extend our comparison with Huang's method using heuristically 



adapted patches introduced in Section 5.2 in order to reduce the computation 



domain to be around the shape contours. To be specific, we sampled along 
the shape contour at an equal geodesic distance d (usually 5 < d < 20), 
and placed a patch centered at each sample point with its scale set to be 
Tp = max (v m , Kj * 

n©( x )(p)) , where r m = pd was the minimum radius so that 
the patches overlapped with one another (1.5 < p < 3, and n = 2). In other 
words, the patches' scale increased proportionally if the two shapes were far 
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Figure 14: Registration using adaptive patches (hand), a) Results of Huang's method 
(mean: 0.29, max: 2.00, variance: 0.22). b) Adapted partitions and results (mean: 0.22, 
max: 2.00, variance: 0.19). g) Regularly distributed patches, c) Regular partitions and 
results (mean: 0.22, max: 2.00, variance: 0.18). 



from each other. Figure 14, Figure [15] and Figure [16| show examples of shape 
registration using adapted patches. In all these examples, the adapted-patch 
method achieved results comparable to the ones using regular partitions and 
better results than Huang's method, while reducing the computational time 
by at least 50 percent. More reductions in computational cost are possible if 
we use finer partitions based on prior knowledge provided, for example, by 
statistical models of specific shapes. 
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Figure 15: Adapted partitions (bunny), a) Results of Huang's method (mean: 0.48, max: 
3.16, variance: 0.39). b) Adapted partitions and results (mean: 0.40, max: 2.83, variance: 
0.30). g) Regular partitions and results (mean: 0.40, max: 2.83, variance: 0.29). 



7.3. Cell division 

Finally, we demonstrate an application of our method on analyzing cell 
deformation. Tracking cell contour is of great interest in biomedical image 
analysis. Many existing methods are based on active contours and level-set 
representations. Although these methods produce good results in segmenting 
cell contours, few of them provide correspondences between deforming con- 
tours. Our method can be combined with level-set segmentation algorithms 



to further establish these correspondences. Figure [17] shows a sequence of 
cell mitosis. The top row shows the segmentation results obtained using the 
level-set segmentation method by Chunming Li et al. [33]. The second row 

51 




Figure 16: Adapted partitions (corpus callosum). a) Results of Huang's method (mean: 
0.36, max: 2.0, variance: 0.25). b) Adapted partitions and results (mean: 0.25, max: 1.41, 
variance: 0.19). g) Regular partitions and results (mean: 0.21, max: 1.00, variance: 0.17). 

shows overlaid shape contours from adjacent images. The third row shows 
the correspondence established using our method with adapted partition-of- 
unity. The deformed source contour (green) is well aligned with the target 
(blue). Interestingly, at the moment of cell mitosis, shape topology changed 
but our method still established reasonable correspondence between shapes 
of different topologies, considering that no continuous deformation actually 
exists to change a shape's topology. 
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Figure 17: Cell-mitosis sequence. Row 1: Cell contours using level-set segmentation in [33] . 
Row 2: Overlaid source (blue) and target shapes (red). Row 3: Correspondence established 
using our method. Row 4: Target shape overlaid over the deformed source shape (green). 

8. Conclusions and future work 

A meshless nonrigid shape-registration algorithm was presented. The reg- 
istration functional is a variational extension of the classic chamfer-matching 
energy in which distance transforms provide registration-error gradients, fa- 
cilitating efficient registration. Also, we modeled shape deformation using a 
meshless parametric representation. This model does not rely on a regular 
control-point grid, and can be adapted to arbitrary shapes. Thus, registration 
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can be focused around the shape contours, greatly improving computational 
efficiency. We tested the proposed method by registering a number of syn- 
thetic shapes, and a deforming cell sequence. Future work includes a 3-D 
extension of the method, and the handling of topological changes. 

Despite promising results, our method still encounters problems in reg- 
istering shapes that have large curvatures, and undergo high-degree defor- 
mation, causing local minima in the registration error. We believe that this 
problem can be addressed by adopting global-optimization algorithms such 
as simulated annealing [5]] , or by the inclusion of statistical priors [6] . 
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